Self-trapping of a Fermi super-fluid in a double-well potential in the BEC-unitarity 

crossover 
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We derive a generalized Gross-Pitaevskii density-functional equation appropriate to study the 
Bose-Einstein condensate (BEC) of dimers formed of singlet spin-half Fermi pairs in the BEC- 
unitarity crossover while the dimer-dimer scattering length a changes from to oo. Using an effective 
one-dimensional form of this equation, we study the phenomenon of dynamical self-trapping of a 
cigar-shaped Fermi super-fluid in the entire BEC-unitarity crossover in a double-well potential. A 
simple two-mode model is constructed to provide analytical insights. We also discuss the consequence 
of our study on the self-trapping of an atomic BEC in a double-well potential. 
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I. INTRODUCTION 



After the experimental realization of Bose-Einstein 
condensate (BEC) and its controlled study under differ- 
ent trapping conditions there have been many inter- 
esting experiments with a cigar-shaped BEC in a quasi 
one-dimensional (ID) trap with a tight transverse con- 
finement Along the axial direction several differ- 
ent types of traps have been employed: harmonic 
double- well 0], periodic optical-lattice [J], and bichro- 
matic optical-lattice [f| traps. Many novel phenomena 
have been predicted and observed in such quasi- ID set- 
ting. Of these, the ones worth mentioning include the 
formation of bright @, gap 0, and dark [§| solitons, 
self-trapping 0, [3] and Josephson oscillation [3, E3, EH ■ 

Macroscopic dynamical self-trapping and Josephson 
oscillation were predicted theoretically [3, E3, El, E3, EH 
and observed experimentally 13, [Tol. Josephson effect was 
observed in super-fluid (SF) 3 He[3 and 4 Hc 0. After 
the experimental observation of BEC in a optical-lattice 
trap [1(| , controlled studies of Josephson oscillation and 
self-trapping in a cigar-shaped BEC seems well under 
control [3J. The studies of such phenomena in a cigar- 
shaped BEC usually employ a double- well potential [||. 
In the simplest case of such a symmetric ID potential 
with the origin of the axial coordinate x set at the trap 
center, under certain initial conditions, when a BEC is 
released with a population imbalance between two sides 
of x = 0, it executes undamped Josephson oscillation on 
both sides of the trap center maintaining a time-averaged 
population imbalance equal to zero. Under different ini- 
tial conditions, the BEC exhibits self-trapping, occupying 
preferably one side of the trap, thus maintaining a def- 
inite non-zero value of time-averaged population imbal- 
ance. The understanding of the transition from Joseph- 
son oscillation to self-trapping and vice versa has been 



*adhikari@ift. unesp.br; URL: www.ift.uncsp.br/users/adhikari 
thpu@ricc.edu 



the topic of many recent investigations [3] . 

A SF Fermi gas in a double-well potential is perhaps 
even more interesting, nevertheless much less studied 
[lc| . (There have been studies of Josephson oscillation 
of a Fermi gas in a OL potential [H[). Such a trapped 
SF Fermi gas gives us the unique opportunity to study 
the Bardecn-Cooper-Schricffcr (BCS) to BEC crossover 
in a two-component Fermi gas under an entirely different 
set-up. The BCS-BEC crossover can be realized by vary- 
ing the attraction between the spin-half fermions forming 
pairs using the Feshbach resonance technique. As the at- 
traction is increased from zero, the simple BCS SF turns 
into a complex Cooper-pair-induccd strongly interacting 
SF and at unitarity (when the Fermi-Fermi scattering 
length cif tends to infinity), it is possible for the Cooper 
pairs to turn spontaneously into Fermi dimers (two-body 
bound state of fermionic atoms) and the BCS SF turns 
into a BEC of dimers. 

After the experimental realization [53] of the BCS- 
BEC crossover in a trapped Fermi SF by varying the 
atomic interaction near a Feshbach resonance, there have 
been renewed interests [2l|, [22|, [23| in the study of a Fermi 
SF at unitarity and beyond in the BEC region where 
we have the BEC of dimers. One can thus recover the 
bosonic behavior in the BEC limit of the crossover (when 
the dimcr-dimer scattering length a tends to zero), while 
expecting new and distinct behavior in the vicinity of the 
unitarity regime. Moreover, on the experimental front it 
is easier to realize a controlled BEC-unitarity crossover 
(BEC side of the BCS-BEC crossover) of the Fermi SF 
than the BCS-unitarity crossover (BCS side of the BCS- 
BEC crossover), as the super-fluid transition temperature 
in the BCS side of the crossover is very low and difficult 
to achieve. 

Here we present a unified Galilean-invariant dynamical 
equation for the study of the BEC-unitarity crossover of 
a cigar-shaped BEC of dimers formed of Fermi atoms. 
In the BEC limit of small dimer-dimer scattering length 
a, the present equation reduces to the usual Gross- 
Pitaevskii (GP) equation [24[ for bosons, and in the uni- 
tarity limit it yields a density-functional (DF) equation 
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|25| for fermions. Hence we call this equation a DF 
GP equation for a Fermi SF valid in the BEC-unitarity 
crossover. For the study of a cigar-shaped Fermi SF along 
the BEC-unitarity crossover, we reduce the present DF 
GP equation to a quasi-lD form and use this reduced 
equation to the study of the self-trapping of Fermi SF 
in a double-well potential. This analytic development 
is presented in Sec. [TTJ This reduced equation also de- 
scribes an atomic BEC with repulsive atomic interaction 
in the BEC-unitarity crossover with different numerical 
value(s) of certain parameter(s), and hence the results of 
the present investigation are also applicable to the self- 
trapping of a repulsive BEC in a double-well potential. 

The numerical simulation with the time-dependent 
quasi-lD equation for a cigar-shaped Fermi SF in a sym- 
metric double-well potential with an initial population 
imbalance between two wells reveals interesting features 
of the Josephson oscillation and self-trapping across the 
BEC-unitarity crossover. In the limit of zero nonlincarity 
one has AC Josephson oscillation. As nonlincarity is in- 
creased by increasing the dimer-dimer scattering length 
or the number of particle, the Josephson oscillation stops 
and self-trapping emerges for a double- well potential with 
appropriate parameters. With further increase of nonlin- 
earity self-trapping is destroyed and the population in the 
two wells executes irregular oscillation. For very large 
nonlincarity, however, the regular Josephson oscillation 
comes back. Nevertheless, for a small number of par- 
ticles, the critical nonlinearity required for one of these 
phenomenon may not be attained and that particular 
phenomenon may not be realized. (The nonlincarity ac- 
tually saturates for a large value of a and hence cannot 
be arbitrarily increased by increasing a as one approaches 
unitarity for a small number of atoms.) These features 
are discussed in detail in Sec. IIVI where we present the 
numerical results. 

In Sec. IIIII we present a simple analytic two-mode 
model to understand the essential features of the numer- 
ical results reported in Sec. IIVI and also point out the 
limitation of the two-mode model. Finally, in Sec. IVl we 
present a brief summary and conclusion of the present 
investigation. 



II. DF GP EQUATION FOR A FERMI SF IN 
THE BEC-UNITARITY CROSSOVER 

At unitarity the following density-functional (DF) 
equation for trapped SF fermions [2(| [27j has produced 
results for energy in clos e ag reement with independent 
Monte-Carlo calculations 1281 



^—V A + U + ^{n) -i — 



2m 



Of 



*(r,t)=0, 



(1) 



where U is the trapping potential, m is the mass of a 
dimcr (twice the atomic mass), ^t(n) = £h 2 n 2 / 3 /m is 
the bulk chemical potential of dimers with density (of 



dimers) n — \^\ 2 , and £ = 2(6n 2 ) 2 ^ 3 (. The normaliza- 
tion condition of the DF wave function is J \*B\ 2 d 3 r = N, 
where N is the number of dimers. At unitarity the only 
length scale is n~ 1//3 , and from dimensional argument the 
chemical potential — and all energies of the trapped SF 
fermions — have the above universal form 12911. 

There have been many theoretical [30L l3l| and ex- 
perimental [32l | investigations which extracted the value 
of the constant £ for fermions, and the most accurate 
value of this constant is given by independent Monte- 
Carlo calculations by two groups [3(|: C ~ 0.44, con- 
sequently, £ sa 13.37 and we shall use this value of £ 
in the present study. For a trapped atomic BEC, the 
energy and chemical potential have the same universal 
form: ~ ^h 2 n 2 / 3 /m [331 ] . where now mass m and den- 
sity n refer to bosons. For fundamental bosonic atoms, 
microscopic numerical calculation based on Jastrow vari- 
ational wave function yielded for the constant £ a slightly 
different value £ = 22.22 [33|. With this modification in 
the value of £, the present investigation could be applied 
to the study of self-trapping of an atomic BEC. 

At unitarity the Fermi pair can stay as a Cooper pair 
or a dimer and transform into each other without transfer 
of energy and Eq. (fTJ can describe both the Cooper pair 
and dimer phases. Here we interpret Eq. (TTJ as the 
equation for dimers. At unitarity the scattering length 
a of two dimers goes to infinity: a — > 00. (Actually, 
at unitarity, the scattering length of two Fermi atoms 
af — > 00. Model studies have indicated that o«a/ (34[. 
Consequently, at unitarity we take a — > 00.) 

Although Eq. (TTJ describes both the dimer SF and 
the Cooper-pair induced BCS SF at unitarity, the bulk 
chemical potential /x(n) appearing in this equation should 
be interpreted differently. For the BCS SF it originates 
from the kinetic energy of Fermi atoms put in differ- 
ent quantum orbitals consistent with the Pauli principle 
discounted for by the negative attractive energy due to 
atomic interaction. For the dimer SF it originates solely 
from the repulsive interaction energy among dimers. In 
the BCS limit as af — > CP wc have the finite nonlinear 
term fj,(n) = 2E F = 2{&TT 2 ) 2 / 3 h 2 n 2 / 3 /m [H in Eq. Q 
originating from the kinetic energy of Fermi atoms with 
negligible contribution from inter-atomic attraction. On 
the other hand, in the BEC limit as a — > + the nonlin- 
ear term for dimers reduces to zero and at unitarity the 
nonlinear term in Eq. (TTJ originates from the saturation 
of repulsive dimer-dimer interaction as a — > 00. 

For the Fermi SF of dimers (and also for an atomic 
BEC) in the BEC-unitarity crossover the following two 
leading terms of the bulk chemical potential of a dilute 
uniform gas can be obtained [351 ] from the expression for 
energy per particle as obtained by Lee, Huang, and Yang 

M 

fjt{n, a) = (4-Kh 2 an/m) (l + a(n 1/3 a) 3/2 + , (2) 

where a = 32/(3-^7?) and n 1//3 a is the dimensionless gas 
parameter. In this expression the scattering length a 
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must be positive (a > 0) corresponding to a repulsive 
interaction. Higher-order terms of expansion @ has also 
been considered [H?}; the lowest order term was derived 
by Lenz [38j | . Considering only the lowest-order term in 
expansion |2[). appropriate in the BEC limit as a, a/ — ► 
+ , the dimers obey the usual GP equation [24[ 
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Considering the second term in the expansion @, in 
the BEC limit, the following modified GP equation for 
dimers can be written following the suggestion of Fab- 
rocini and Polls [391 
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Equation (0J provides an adequate correction to the GP 
equation ([3]) for small a. But as a increases and diverges 
at unitarity, the nonlinear term should saturate to the 
finite universal nonlinear term /J,(n) of Eq. |T]) and not 
diverge like the nonlinear terms of the GP equation Q 
and of the Fabrocini-Polls equation (j4]) . The chemical po- 
tential and energy should not diverge at unitarity, as the 
interaction potential remains finite in this limit, although 
the scattering length a diverges. In the weak-coupling 
GP limit, the scattering length serves as a faithful mea- 
sure of interaction. But as a increases, it ceases to be a 
measure of interaction. For a general scattering length, 
an exact expression of the chemical potential is not avail- 
able. However, a recent quantum Monte Carlo study 
maps out the equation of state in the entire BEC-BCS 
crossover regime f30j . 

Following a recent suggestion [i(J, for the full BEC- 
unitarity crossover we consider the DF GP equation for 
the dimer SF providing a smooth interpolation between 
Eqs. (d|) and Q: 
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l + a<(l + <5)a 3 / 2 |*| 



to (an 1 / 3 ) be continuous at unitarity. The condition for 
continuity yields a small value for <5(~ 0.04). How- 
ever, we shall take 5 = in this study. This will make 
further analytical development easier while maintaining 
the first derivative of fi(n,a) with respect to (an 1 / 3 ) ap- 
proximately continuous at unitarity. A set of equations, 
similar to Eqs. |5|) and (JB]), for fundamental bosons, and 
not for composite dimers, produced results for energy 
[2(| of a trapped condensate in agreement with Monte- 
Carlo calculations |4j|. A similar equation for the BCS- 
unitarity crossover produced results for energy [26|, [4l[ of 
a trapped BCS SF in agreement with Monte-Carlo cal- 
culations 44[ . Different parametrization of the chemical 
potential in the BEC-BCS crossover have been proposed 
in the literature [45| . We do not expect our results in this 
work will be sensitive to which specific form for /x(n, a) 
we choose to use here. Furthermore, it has been shown 
that the DF GP equation ([5]) is equivalent to the quan- 
tum hydrodynamic equations for dimers p6l . [45| 
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where s is a phase, and v is the velocity. 

For a cigar-shaped SF, where the transverse trapping 
is very strong, the interesting dynamics is confined in 
the axial direction and in the transverse direction the 
system is confined in its ground state. In such a quasi-lD 
geometry, the axial and transverse coordinates decouple 
and it is useful to write an effective ID equation for the 
dynamics of a cigar-shaped SF and we perform the same 
in the following. For the cigar-shaped double-well trap, 

U(r) = V(x) +mu; 2 (x 2 + X 2 z 2 + X 2 y 2 )/2, 
V(x) = Aha; exp(—Kmuix 2 /H) , 

where A ^> 1, and A and n are two di- 
mensionless parameters characterizing the strength 
and width of the barrier, respectively, it is ap- 
propriate to take *(r, t) = ip(x,t)(j)(y)(j)(z) with 



1 + a<5a 3 / 2 |*| + 7 (1 + <5)a 5 / 2 |*| 5 / 3 ) ^ ^ = ["^VM 1 / 4 exp[-ma;Aj/ 2 /(2^)] 



where n = |*| 2 and 8 and 7 are yet unknown constants 
satisfying 7 = Ana/^. (These equations arc also valid for 
an atomic BEC with m and a representing atomic mass 
and scattering length, respectively, and with £ = 22.22 
[HI.) By construction, Eq. © yields the limit (J2J for 
small a; it also has the correct behavior at unitarity. Us- 
ing a similar expression for /x(n, a) in the BCS-unitarity 
crossover [4l[ the constant S was calculated [42[ by re- 
quiring that the first derivative of /x(n, a) with respect 



representing 

the harmonic-oscillator ground state in the transverse di- 
rection and ij}(x,t) representing the essential dynamics 
in the x direction. The potential V(x) together with the 
harmonic trap mu> 2 x 2 jl simulate a double- well in the ax- 
ial x direction. Multiplying Eqs. (fTJ) and ([4]) by (j>(y)<j)(z) 
and integrating over y and z we get the following ID 
equations [27l. |46| 
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where n = \ip(x, £)| 2 , U{x) = Ae~ KX ~ + x 2 /2 represents 
the double- well potential and n = \ip\ 2 and we use har- 
monic oscillator dimensionless units h = m = lu = 
1. All lengths are now expressed in oscillator unit 
\Jh/(rnuj) and time in w -1 , and ip is normalized as 
! co oo dx\^x,t)\ 2 =N. 

A simple DF GP equation interpolating between 
Eqs. © and © is 

Id 2 d 



2gx2+U (x)+„(n,a)-i gt 



fi{n, a) = 2aXn 
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where n = \ip\ 2 and = 8aA 5 / 6 7r 1 / 6 /(30- Equation © 
reduces in the BEC a — > + limit to Eq. © and in the 
unitarity a — > oo limit to Eq. ([7]). We shall use Eq. © for 
the description of self-trapping and Josephson oscillation 
in the BEC-unitarity crossover of fermions. 



III. TWO-MODE MODEL OF THE FERMI SF 
A. Two-Mode Model 

Before we present the full numerical results, it is in- 
structive to consider the so-called two-mode model @ 
which is widely used in the study of BEC in a double- well 
potential, and more recently, has been used in the inves- 
tigation of Fermi SF across a weak link [53] . Due to its 
simplicity, the two-mode model can provide many useful 
insights. Here we construct the corresponding two-mode 
model of the Fermi SF based on Eq. © . 

To this end, we decompose i(>(x,t) as 



(11) 



where the spatial mode functions <p\ i2 {x) & r e assumed to 
be real, satisfy the orthonormal condition 



4>i(x)(j)j(x)dx = Sij , 



and are localized in each of the two wells, respectively. 
^1,2 (i) are in general complex and satisfy the condition 
\^'iMt)\ 2 = N 1;2 (t) so that 

\i> 1 (t)\ 2 + \tp 2 (t)\ 2 = N 1 (t)+N 2 (t) = N. 

Inserting the decomposition (fTTj) into Eq. © , integrating 
out the spatial degrees of freedom, we obtain the follow- 
ing equations of motion for ip\, 2 (t): 

#i = Exipx +£i(|Vi|)V'i-£V>2, (12) 
# 2 = 3 > Vi+£a(hM)lfe-£V'x, (13) 
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(15) 



K = - 



+ U(x) Uj(i), (16) 



with rii = \ipi4>i\ 2 - Here we have neglected integrals in- 
volving spatial overlaps of 4>\(x) and 4> 2 (x). 

For simplicity, we assume a symmetric double well with 
U(x) = U(—x) so that 4>i(x) = (f) 2 (—x) and consequently, 
we have Ex = E 2 and £i(|V>|) = ^(M) = £(|Vj )- Le t us 
write the waves il>x,2 in terms of its amplitude yN\~2 and 
phase {6\, 2 ) 



i>x, 2 = 1,2 

and define a pair of conjugate variables: 

S=(N 1 -N 2 )/N, = 02-01. 

Here the variable S denotes the population imbalance 
between the two wells and 9 is the phase difference. After 
some straight-forward algebra, the following equations of 
motion for S and 9 can be derived from Eqs. (fT2)) and 
(p~3|) by equating the real and imaginary parts of both 
sides: 



S = -2K,y/l~S 2 sin0, 
= \E(y/Wx) - £(y/Ni) 



2K- 



S 
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These are the two-mode equations for the Fermi SF. 
Note that the two-mode equations describing weakly- 
interacting bosons in double-well potential 0] are recov- 
ered if we take /i(n, a) = 2a\n and, correspondingly, 
£{y/Wi) = 2aXN l J\<j>,\ 4 dx. 

The two-mode equations can be cast into the canonical 
form 

£, = _dH - OH 

89 ' as" 

with the classical Hamiltonian defined as 



H = 



£{^/n[) - £ (y/N^) dS -2/Cv 7 !- S 2 cos0. 



(17) 

By studying the properties of this Hamiltonian, we can 
tell whether the system should exhibit self-trapping or 
Josephson oscillation. 



B. Fermi SF at Unitarity 

As a concrete example, let us consider the Fermi SF at 
unitarity where a — > oo and from Eq. (|10[) we find 
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FIG. 1: (Color online) Equi-energy phase contour plot of the 
unitary Fermi SF. (a) Contours of different energies for A = 8. 
(b) Contours of E = 10 for different values of A. Energy is in 
units of 2/C. 



It follows that 
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with U = (3£/5)(A/tt) 2 / 3 J|0 4 | lo/3 d2:. The classical 
Hamiltonian takes the form: 
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(i + s) 5/3 + (i-sf /3 -^i-s 2 cose 



where A = {N/2) 2 / 3 U/ (2/C) measures the ratio of the 
strength of the nonlinearity (N/2) 2 ' 3 U and the tunneling 
energy 2/C. 

As can be seen from Fig. [TJ if we draw equi-energy 
contours of H in the phase space of (S,9), we can see 
two types of contours: those that form closed loop and 
those do not. The division of these two types of contours 
occurs at the critical energy (we take 2/C as the units for 
energy) 
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If the system has energy E < £7 C ritj its dynamics will 
follow the closed contours, and both 5 and 9 will oscil- 
late in time. In particular, the population imbalance 5 
oscillates around and the time averaged population im- 
balance vanishes. This corresponds to the AC Josephson 
regime. On the other hand, if the system has energy 
E > -Ecrit, it will follow the open contours where 6 will 
grow indefinitely and 5 will oscillate around a non-zero 
value and will never cross the 5 = line. This corre- 
sponds to the self-trapping regime. 

Given the initial values 5(0) and 9(0), the system 
moves on a contour of constant energy given by 
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3A r 
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(l + 5(0)) 5/3 + (l-5(0)) 5/3 
v/l-5(0) 2 cos6>(0). 
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may be recast into the form 

1+ v/1- 5(0) 2 cos 6>(0) 



A > A r = - 



3(l + 5(0)) 5 / 3 + (l-5(0)) 5 / 3 



(18) 



In other words, within the two-mode model, the ratio 
of the nonlinear strength and the tunneling energy, A, 
determines whether the system should be self-trapped or 
not. 

To determine the values of these quantities, we need 
to choose properly the spatial mode functions 4>i.2{x)- A 
reasonable choice is given by 0, l48j 



0i, 2 (a;) 



<j>+(x) ± </>_(x) 
V2 



where the normalized wave functions </>± (x) are the 
lowest-energy symmetric and antisymmetric stationary 
solutions to the time-independent DF equation: 



1 d 2 
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with chemical potential /i±. 

In Figs. [2](a), (b), and (c) we illustrate the nonlinear 
strength (N/2) 2/3 U, tunneling energy 2/C and their ratio 
A as functions of NX, respectively. As NX increases, the 
strength of the (repulsive) nonlinearity increases. As a 
result, the <j>\ and <f>2 become widened and enjoy more 
overlap. This leads to an increased tunneling energy. 
However, the ratio of the nonlinear strength and the tun- 
neling energy A does not have a monotonic behavior as 
iVA is increased. As shown in Fig. [2jc), A initially in- 
creases for small NX, reaches a peak and then decreases. 
In a recent study, Salasnich et al. [13] used the local- 
density approximation on top of quantum Monte-Carlo 
data of Ref. [3(| to explore the phase diagrams and find 
regimes of Josephson tunneling and of dynamical self- 
trapping of a 3D Fermi supcrfluid. In the two-mode ap- 
proach reported in Ref. [47| , a constant tunneling energy 
is arbitrarily chosen for the whole crossover regime. This 
is an inappropriate oversimplification. In Fig. [UJd), we 
show the critical value A c as a function of the initial 
population imbalance 5(0). One can see that as 5(0) 
increases, A c decreases rapidly. 

The fact that A is bounded from above even though 
the interaction strength can increase without bound has 
important consequences. For example, for certain ini- 
tial conditions, self-trapping may only occur within an 
intermediate range of nonlinearity. Both too small and 
too large a nonlinearity will destroy self trapping. This 
statement is also true away from the unitarity, even in 
the BEC limit. Furthermore, for a sufficiently small 5(0), 
A may never exceed the corresponding A c . When this is 
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FIG. 2: (a) Strength of nonlinearity (N/2) 3/2 U, (b) tunneling 
energy 2/C, and (c) their ratio A (c) as functions of iVA at 
unitarity. (d) The critical A c as function of 5(0) for 8(0) — 
[see Eq. (fit))]. 



the case, the system will always stay in the Joscphson 
oscillation regime. For example, according to Fig. f^d), 
A c w 300 for 5(0) = 0.1. Fig. ^c) shows that the sys- 
tem can therefore never reach the self-trapping regime if 
5(0) equals 0.1 or smaller. This is consistent with our 
numerical results. 

We want to remark that even though results obtained 
from the simple two-mode model may provide significant 
qualitative insights, they are not expected to be accurate 
quantitatively. Particularly for large nonlinearity, predic- 
tions from the two-mode model can deviate greatly from 
the numerical results [48|, [4^ . The error mainly occurs 
in estimating the tunneling energy IC. The two-mode 
equations (|12j) and (|13p are obtained by neglecting many 
terms involving overlap integrals of the mode functions 
4>\ and 4>2 , and hence in general greatly underestimates 
the tunneling energy, particularly for large nonlinearity 
when overlap between <pi and <p2 can be significant. Fur- 
thermore, the two-mode approximation itself becomes 
questionable for large nonlinearity. When there is ex- 
change of atoms between the two wells, the mode func- 
tions will change accordingly due to the modification of 
the nonlinear mean-field. Indeed, in our numerical calcu- 
lations to be presented below, we observe that the spatial 
wave function of the system changes in time. In certain 
regimes, this change is significantly enough to invalidate 
the two mode model. 



IV. NUMERICAL RESULTS 

In this section we present an account of the numeri- 
cal study of self-trapping and Josephson oscillation in a 
double-well potential by solving the full quasi-lD DF GP 
equation ([9]) valid for a cigar-shaped SF. The double- well 
potential is taken as 



U(x) = x 2 /2 + Ae~ 



(19) 



We shall take the parameters A and k of this double well 
similar to the ones employed in Ref. (4fj| in a study of 
self-trapping with dipolar bosonic atoms. 

To create an initial state with desired population im- 
balance for a given set of parameters NX and a, we search 
for the ground state of an asymmetric well comprised of 
an off-centered harmonic potential and the Gaussian bar- 
rier potential 



U'(x) = (x - x f/2 + Ae 



(20) 



The ground state of this asymmetric well is obtained by 
solving the time- independent version of the DF GP equa- 
tion ([9]) using the imaginary time evolution method. The 
parameter xq in (|20p is chosen so that the population im- 
balance 



S(t) = (N 1 (t)-N 2 (t))/N. 



(21) 



has a fixed pre-determined initial value 5(0). Here Ni(t) 
and N2(t) are the number of dimers in the first and the 
second well of the double-well potential. Experimentally, 
this is indeed the method used to generate the initial 
population imbalance [|[ . We have also considered other 
forms of initial wave functions and found that the final 
results are qualitatively insensitive to the specific forms 
provided the initial population imbalance 5(0) is kept 
fixed at a small value. However, at a quantitative level 
the results could be sensitive to the form of the initial 
wave function. The sensitivity of the result to the initial 
wave form increases as 5(0) is increased. Moreover, the 
results are quite sensitive to the initial 5(0) employed. 

Once the initial wave function is chosen, Equation ^ 
is solved numerically after discretization with the Crank- 
Nicolson scheme (50, [EH employing space and time steps 
0.025 and 0.0002, respectively, using real-time propa- 
gation with the FORTRAN programs provided in Ref. 
[501 ] . The results are also independently confirmed using 
a MATLAB code based on the split fast Fourier trans- 
form method. 

Now we present results of dynamical evolution of a 
Fermi SF, where we take £ = 13.37 in Eqs. © and (flO|> . 
The numerical study of self-trapping and Josephson os- 
cillation with the Fermi SF of dimers along the BEC- 
unitarity crossover reveals interesting features. To start 
the investigation of self-trapping we fix the trap param- 
eters A and k in Eqs. (fT^|) and (|2T)|) at nontrivial values 
(a not too small value of A and a not too large k) , which 
permit smooth and free Joscphson oscillation in the BEC 
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FIG. 3: (Color online) Population imbalance S(t) vs. t dy- 
namics with (a) AA = 100, 5(0) = 0.2, k = 10, A = 8, (b) 
iVA = 100, 5(0) = 0.2, k = 10, A = 16, (c) 7VA = 1000, 5(0) = 
0.3, k = 10, A = 16, (d) iVA = 100, 5(0) = 0.3, k = 8, A = 12, 
for dimer-dimer scattering length a varying over the BEC- 
unitarity crossover. 



limit (a = 0). Note that a very small value of A and a 
very large value of k tend to reduce the double well l|19p 
to a single well where there cannot be any self trapping 
and Josephson oscillation should appear for all values of 
a and N. In order to have self trapping, A cannot be 
too small and k cannot be too large. This is illustrated 
in Fig. [1(a) for A = 8,k = 10, NX = 100,5(0) = 0.2 
and for a = 0.1,0.01, and 0.001 where we plot S(t) vs. 
t. There is no self trapping for a very small value of 
S(0)(= 0.1). The quantity S(t) is experimentally mea- 
surable and S(t) vs. t dynamics provides information 
about self trapping and Josephson oscillation. From Fig. 
[3] (a) we find that there is Josephson oscillation for all val- 
ues of a and there is no sign of self trapping. (A nonzero 
time average (S(t)} ensures self trapping.) But a com- 
pletely new scenario emerges as A is increased to 16 from 
8, as can be seen from Fig. [3] (b) where we show the S(t) 
vs. t dynamics for a = 0.01,0.1,0.2 and 0.5. The plots 
for a = 0.01 and 0.5 of Fig. [3] (b) are quite similar to 
the plots for a = 0.001 and 0.1 of Fig. [3] (a) illustrating 
regular (periodic) Josephson oscillation with no sign of 
self trapping. But, for intermediate values 0.1 and 0.2 of 
a, self trapping and irregular (non-periodic) oscillation 
can be seen in Fig. [3J (b). In Figs. [3J (c) and (d) we 
illustrate two more cases of S(t) vs. t dynamics with a 
different value of 5(0) (= 0.3) and different NX and trap 
parameters, respectively, where one can clearly find self 
trapping. 

In the following, we discuss in detail the results for 
three initial population imbalance 5(0) = 0.1, 0.2 and 
0.3, which are representative for a general case. 

Population imbalance 5(0) = 0.1 : For this relatively 
small initial population imbalance, we found that for any 
values of NX and a, the system is always in the Josephson 
regime: the population imbalance S(t) oscillates sinu- 
soidally between 5(0) = -0.1 and 5(0) = 0.1. The sys- 
tem never exhibits self-trapping. The frequency of oscil- 
lation increases as the strength of nonlincarity increases. 
Note that the strength of nonlincarity is increased by in- 
creasing either A^A or a. However, the nonlinear interac- 
tion among dimers saturates as scattering length a — ► oo 
at unitarity, it increases indefinitely with A^A. This re- 
sult is consistent with our previous discussion of the two- 
mode model: For a sufficiently small 5(0), the required 
critical value of A c for self-trapping cannot be achieved 
by increasing the strength of the nonlincarity and the 
system stays in the Josephson regime for all values of 
A^A and a. 

Population imbalance 5(0) = 0.2 : The results for the 
S(t) vs. t dynamics for this initial population imbalance 
is illustrated in Figs. [3] (a) and (b) for fixed A^A = 100 and 
various values of a and two different traps as described 
earlier. In Fig. [3] (b), for a small scattering length of 
a = 0.01 (solid line), Josephson oscillation is observed. 
When a is increased to 0.1 (dashed line), self-trapping is 
clearly seen — S(t) does not deviate from 5(0) too much 
and never crosses zero. With further increase of a to a 
slightly larger value of 0.2 (dotted line), self-trapping is 
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increasing nonlinear interaction strength 

FIG. 4: Different dynamical regimes of Fermi super-fluid in 
the BEC-unitary crossover. 



destroyed and S(t) exhibits irregular oscillations around 
zero. Accompanied with this irregular population oscilla- 
tion, the density profile \ijj(x,t)\ 2 also develops complex 
and irregular structures. Remarkably, upon further in- 
crease of a, as the dot-dashed curve for a = 0.5 shows, 
regular oscillation returns and the population dynamics 
once again exhibits sinusoidal Josephson oscillations. 

Population imbalance 5(0) = 0.3 : Finally, let us dis- 
cuss this relatively large initial population imbalance. If 
we use NX = 100 as we did above for 5(0) = 0.2, when 
a is increased from zero to oo, the system sequentially 
makes transitions from Josephson, to self-trapping and 
finally to irregular oscillation regimes. The Josephson 
oscillation is never recovered for NX = 100 for very large 
values of scattering length a (results not shown here) . In 
Fig. [3] (c) we plot the results for 5(t) vs. t dynamics for 
NX = 1000. In this case, in addition to the three regimes 
just mentioned, for a sufficiently large a(= 10), Joseph- 
son oscillation is restored, just as in the case of 5(0) = 0.2 
and iVA = 100 discussed above. In Fig. Fig. [3] (d) we 
show another example of S(t) vs. t dynamics for a differ- 
ent trap, which is quite similar to that in Fig.[3j(c). We 
also did some calculation with larger 5(0) where a sim- 
ilar panorama emerges and we do not report the details 
here. 

To summarize the general characteristics of the pop- 
ulation dynamics, we find that for any given initial pop- 
ulation imbalance and for either sufficiently small or suf- 
ficiently large nonlinear interaction strength, the system 
is in the Josephson oscillation regime. For intermedi- 
ate interaction strength, the system can make transi- 
tion to self-trappping and irregular oscillation regimes 
as schematically shown in Fig. 5] The critical interac- 
tion strength at which the system makes the transition 
to self-trapping is sensitive to the initial population im- 
balance and increases sharply as 5(0) increases. (It is 
also sensitive to the parameters for the Gaussian barrier 
that creates the double- well potential.) It is possible that 
for a sufficiently small 5(0), the system always stays in 
the Josephson regime. The restoration of the Josephson 
oscillation at large interaction strength may seem sur- 
prising at first sight. However, one can understand it in 
the following intuitive way. For a sufficiently large inter- 
action strength, the chemical potential is large and the 
effect of the Gaussian barrier becomes relatively unim- 
portant. The wave functions on opposite sides of the 
barrier have sufficient overlap and hence the cloud tun- 
nel back and forth without difficulty. 

The appearance of self-trapping is best illustrated 
through a study of the time-averaged population imbal- 
ance (5(f)) vs. nonlincarity aNX and we do that next. In 
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FIG. 5: (Color online) Time-averaged population imbalance 
(S(t)) vs. nonlinearity aNX for different NX for trap param- 
eters Al = 16, k = 10 and (a) 5(0) = 0.3 and (b) 5(0) = 0.6. 
Different curves are generated by varying the scattering length 
a across the BEC-unitarity crossover for S(0) = 0.3 and 0.6, 
respectively, based on initial wave form (1201) . 



Fig. [5] (a) we plot (5(f)) vs. aNX by varying the scatter- 
ing length from to oo for a fixed NX with trap parame- 
ters A = 16 and k = 10. The initial population imbalance 
is chosen as 5(0) = 0.3. In Fig. [5] (a), with the increase 
of aNX, self-trapping appears for aNX slightly greater 
than unity. With further increase of aNX, self-trapping 
increases with an increase of (S(t)). For NX = 10, self 
trapping never disappears and continues even at unitar- 
ity However, beyond aNX « 5 self-trapping decreases 
with the increase of aNX for NX = 20, 100 and 1000. For 
larger NX, (5(t)) eventually goes to zero as the nonlinear 
repulsion becomes too large to maintain all dimers in a 
single trap, except for NX = 10. In Fig. O (b) we ex- 
hibit a similar plot for 5(0) = 0.6 with trap parameters 
A = 16 and re = 10. 

Finally, to check the validity of the quasi-lD approx- 
imation, we performed full 3D numerical simulations 
based on Eqs. ([5]) and The quasi-lD approximation 
should be valid when \i -C Xhoj. We have chosen different 
sets of parameters, some of which satisfy and the rest vi- 
olate the quasi-lD condition. For parameters such that 
the quasi-lD condition is satisfied, we indeed find that 
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our 3D numerical results are nearly identical to the ID 
results presented here. For parameters that the quasi-lD 
condition is violated, the 3D results show deviations from 
the ID ones. However, the qualitative features (i.e., the 
dependence of different dynamical regimes on the initial 
population imbalance and the strength of nonlincarity) 
presented in Figs. [3] and O remain valid. Specifically, we 
verified that the results reported in Figs. [3] remain es- 
sentially valid in the 3D model. This assures us of the 
reliability of the present study in ID. 

V. SUMMARY AND CONCLUSION 

To summarize, we have studied the dynamical proper- 
ties of a Fermi SF confined in a double-well potential in 
the BEC-unitary crossover regime. To this purpose, we 
have developed a nonlinear Schodinger equation valid in 
the whole regime based on a density functional approach 
and on the equations of state from quantum Monte Carlo 
calculations. This equation is equivalent to the hydro- 
dynamic equations with the quantum pressure term in- 
cluded. In the BEC side of the crossover, it describes ac- 
curately the equilibrium and low-energy dynamical prop- 
erties of the Fermi SF. In particular, Josephson effect 
has been investigated using this method [l8| and the re- 
sults have been shown to agree with those obtained from 
the microscopic approach by solving the Bogoliubov-de 
Gennes equations [52[. Compared with the latter, the 
great advantage of the current approach is its mathe- 
matical simplicity. 

We have identified three dynamical regimes of the sys- 
tem: the Josephson regime, the self-trapping regime and 
the irregular oscillation regime. For a given initial pop- 
ulation imbalance, these regimes are accessed according 
to the strength of nonlincarity as schematically shown 
in Fig. 3J The Josephson regime is always reached 
at either sufficiently small or sufficiently large interac- 
tion strength. For a small initial population imbalance, 
Josephson regime may be the only regime that the system 
can have access to. Note that the strength of nonlincar- 
ity can be increased by either increasing the number of 
dimers N or the scattering length a. However, it satu- 
rates as a tends to infinity while no saturation occurs for 
large N. 



The quasi-lD model — Eqs. ([9]) and (JTOj) — presented 
and used in the study of dynamical evolution of a Fermi 
SF in the BEC-unitarity crossover in this paper is also 
valid for an atomic BEC with a slightly modified value 
for the parameter £. Hence the present results for self- 
trapping of a Fermi SF in a double-well potential are also 
applicable for a repulsive atomic BEC when the atomic 
scattering length varies from to oo. However, in this 
case there could be practical difficulty with three-body 
loss in the experimental realization of the system for large 
scattering length. 

We have also developed a simple analytical two-mode 
model, analogous to the much studied system of a BEC 
in a double-well potential. We show that the properties 
of the system can be described by a classical Hamilto- 
nian with population imbalance and relative phase as a 
pair of conjugate variables. The great advantage of the 
two-mode model is its simplicity which makes analytical 
studies possible. The key parameters that characterize 
the two-mode model are the strength of nonlinearity and 
the tunneling energy. We calculated these parameters 
using the spatial mode function obtained by numerically 
solving the full time-independent nonlinear Schrodiongcr 
equation. From this calculation we show that the ratio 
of the interaction strength and the tunneling rate can- 
not increase indefinitely when the interaction strength in- 
creases. This explains the numerical observation that for 
sufficiently small initial population imbalance, the sys- 
tem may always stay in the Josephson regime. However, 
care must be taken when making quantitative compar- 
isons with numerical results. In particular, for strong 
nonlincarity, the two-mode model can be even qualitative 
incorrect. For example, this model predicts the existence 
of the Josephson and the self-trapping regime, but not 
the irregular oscillation regime found in the numerical 
calculation, which occurs at relatively large nonlinearity 
and lies in the regime where the two-mode model is no 
longer valid. 
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